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Employing a general variational method and perturbation theory, we derived explicit solutions for 
the description of one-dimensional two species Bose-Einstein condensates confined by a harmonic 
trap potential in an optical lattice. We consider the system of two coupled Gross-Pitaevkii equations 
(GPE) and derive explicit expressions for the chemical potentials and wavefunctions in terms of the 
atom-atom interaction parameters and laser intensity. We have compared our results with the 
numerical solutions of the GPE and performed a quantitative analysis for the both considered 
methods. We underline the importance of the obtained explicit solutions to characterize the density 
profile or degree of miscibility of the two components. 
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I. INTRODUCTION 

Multiple Bose-Einstein condensates (BEC) of different 
atomic species have been realized in the last years. 
Mixture of alkali atoms of 87 Rb in two different hy- 
perfine internal spin states,® atoms 23 Na with a su¬ 
perposition of spinor condensates J- 2 combination of 41 K 
- 87 Rb|® s7 Rb- 85 Rb,® 87 Rb- 133 Csj® and gases of rare 
atomic species 168 Yb- 174 Ybj®have been employed to pro¬ 
duce two species BEC. These quantum degenerate mix¬ 
tures allow to study several intrigui ng phenomena as 
the dynamics of the superfluid system} 2 ® the production 
of heteronuclear polar molecules p the miscibility or im- 
miscibility of the two quantum fluids, 8 among other ef¬ 
fects. Also, two-sp ecies BEC loaded in a optical lattice 
have been explored!® 10 A similar system but of Fermi- 
Bose quantum gas mixture in a 3-dimensional optical lat¬ 
tice was implemented to study the interfering paths of 
the bosonic wave function scattered by the presence of 
fermionic atoms. 41 These results have led to an intense 
theoretical and mathematical studies on the properties 
of the two-coupled Gross-Pitaevkii equations. 

The basis of this research lies on the knowledge of the 
dependence of the chemical potentials as functions of the 
interparticle interactions and the spatial density proba- 
bilityP 

A fascinating experimental realization to study the 
one dimensional (ID) transport properties of ultracold 
fermionic and bosonic atoms in a periodic potential have 
been reported in Reffl2l 

From the theoretical point of view there are several 
studies for the description of two species Bose conden¬ 


sates. Typically, numerical approaches or Thomas-Fermi 
approximation are employed to calculate the chemical po¬ 
tential and the ground state wave functions. 4 ® In Reffl4[ 
it is analyzed the mixture of ID two interacting conden¬ 
sates modeled by the Bose-Hurbbard Hamiltonian and by 
using the quantum Monte Carlo numerical simulations. 
Theoretical analysis of the ID two component BEC prob¬ 
lem becomes an important reservoir to mimic different 
physical effects of the Condensed Matter Physics (see for 
example Refs|15HT7j), including the magnetic properties 
of the bosonic mixtures with tunable interspecies inter¬ 
actions!- 4 ® Also, as it will be shown below, we can take 
advantage of analytical results for the study of quantum 
effects and predictions for cold atoms researches. 

Assuming a “c igar-shaped” type for the Bose-Einstein 
condensates 19 ^ of a gas composed by two kind of bosons 
loaded in an optical lattice, we can consider the following 
system of ID GP equations: 
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Here, u>i > 0 denotes the harmonic trap frequencies where 
for simplicity we consider the same for both condensates, 
i.e., uq = w 2 = w, nii > 0 , and Vi are, respectively, the 
mass and chemical potential for the specie i (i = 1 and 
2 ), Vl > 0 and d > 0 the intensity and laser wavelength, 
A i takes into account the self-interaction term for the ith 
specie, and A 3 , the interaction between unlike particles 
of the species 1 and 2. In this system, the complex func¬ 
tion <f>, (a:) is knowiP-^ as the macroscopic wavefunction 
or order parameter of the ith component and is defined 
as the expectation value of the corresponding field oper¬ 
ator, namely $i(x) = ( ( P i (x)). The functions satisfy 
the normalization conditions 
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Therefore, the partial Frechet derivatives of E are 

OrE = E'^) + X 3 \M0\ 2 M0 , ( 11 ) 

8 2 E = E'M + X 3 \M0\ 2 M0 ■ ( 12 ) 


The minimum of the energy E{pp\,ip 2 ) under the restric¬ 
tions L |tM £)| 2 = Ni satisfies the Lagrange conditions 

for some constants pi/2 (z = 1 , 2 ), 


/ \§ i {x)\ 2 dx = N i , * = 1,2 , 
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where N t denotes the number of atoms of the zth specie. 

It is worth to notice that in some situations, as in the 
case of spinor condensates, where one produces confine¬ 
ment of an atomic cloud of an element in different spin 
statesp^ll the condition ([5| must be substituted by 

f \$>i(x)\ 2 dx + f |$2 (*)1 2 dx = N , N = Ni + N 2 . 

Jr Jr 

We can rewrite the system (|T|) in its dimensionless 
form, by considering, for instance, l = \Jh/ (toiw), x = 
l£, and $i(a;) = ^(O/v^ * = 1 , 2 , in which case we have 


4 )* + [£ 7 - p] V = 0 , 

where Cq and £/ are respectively the operators 
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Notice that <[13j) coincides with (| 6 j). 

In previous workspHUD we have presented different 
methods to express the chemical potential p and the or¬ 
der parameter ip(£) as function of the interaction param¬ 
eter A for the ID Gross-Pitaevkii equation. In the present 
paper, we adapt two of these methods (the generalized 
variational approactP 5 and perturbation theory) for the 
system ([b]), by considering the vector chemical potential 
p as function of the atom-atom interaction strength of 
each component Ai, A 2 and the interaction between both 
species, A 3 . 

The paper is organized as follows: in Section II we 
present the mathematical framework of the variational 
problem formulation, which characterizes the condensate 
as ground state solution for the system § , as well as 
its equivalent integral representation. We also report an 
exact representation of p( Ai,A 2 ,A 3 ) over which is based 
our variational approach described in Section III. In Sec¬ 
tion IV we develop the perturbation method valid for two 
coupled GP equations. Section V is devoted to present 
the results of these two approaches comparing with the 
exact numerical solution of the system <§■ Also, final 
conclusions are delivered showing the range of validity 
of both considered methods, with respect to parameter 
values employed for the description of two-species Bose- 
Einstein condensate in an optical lattice. 
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II. GENERAL MATHEMATICAL 
FRAMEWORK 


Here, a 2 = mi/m 2 , Aj = A i/lhui, (i = 1,2,3), Vo = 
Vl/Sw, a = 2ht/d and pj = Vj/hx (j = 1,2). For the 
system Q, the energy functional can be cast as 

E^i,^) = Ei(ipi) + E 2 (ip 2 ) 

+Y [ \M0\ 2 \M0\ 2 dt , ( 10 ) 
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In this section we establish the functional framework 
for the mathematical analysis of existence, regularity and 
stability of ground state solutions for the system §. 
There is a great number of mathematical work on these 
questions, some of them mentioned in the references be¬ 
low. The eingenvalue problem ([ 6 ]) has an intrinsic math¬ 
ematical interest, but the ground state solutions (i.e., 


















3 


standing wave solutions of minimal energy) play impor¬ 
tant role for condensates. By standing wave we mean 
solution of the evolution equation 

flip 

i— = [4,+ £,]*, (14) 

of the form 
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A. Existence of ground states and their stability 


The existence of a minimal energy solution is a conse¬ 
quence of the Gagliardo-Nirenberg inequality (see The¬ 
orem 1.3.7 in Ref ETl) . which in ID allows us to show 
that the energy functional E is bounded by bellow on 
the manifold E, for all values of A * Gl,i = 1, 2, 3. With 
arguments of convexity, we can show that the (real) so¬ 
lution of minimal energy is unique provided that Ai, A 2 
and A 3 are positive. Moreover, since the system ([ 6 ]) has 
the properties of conservation of energy and mass (i.e., 
th e num ber of particles), we can prove the orbital stabil- 
it yHHlH of ground states. 

On the other hand, the space T~L = L 2 (M) x L 2 (R) is a 
Hilbert space if one considers the usual inner product 


We consider the following minimization problem 

E min (A) = min{E(\l>); \l> £ E}, (15) 

where A = (Ai, A 2 ,A 3 ), 'E = (^ 1 ,^ 2 ), 
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and S = V x V, where 

V = {$ G H\R) ; f [|^(0| 2 + £ 2 |V>(0| 2 ] d£ < cx)} 

and 7? 1 (M) is the standard Sobolev space. 

Although the solutions of Eq. ([6]) are in general com¬ 
plex valued functions, we can restrict ourselves to just 
the real valued ones. This can easily be justified because 
any solution of this system satisfies the following inequal¬ 
ity!^ there exist 0 < S < 1 and C(5) >0 such that 

l*(0| 2 + l*Wl<^)exp(-^ 2 ), (16) 
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the exponential decay (16) and a simple calculus gives 
TpiRipu - i>xR^'u = ^ l ^ 17 ! 2 = 0 ’ 


Therefore, r ipiR = pijju for some real constant /3 7 ^ 0. 
The same holds for second component of \E, which gives 
us tl> 2 R = 7 ^ 2 / for some constant 7 . Hence, the function 


y /l + P 2 iPir{Q 

Vi +1 2 ^2r{C) 


is a 


u(0 = 

real solution of § and ( [It] ) is given by 

V’i r(0 
i>2 r(0 


(18) 


*(0 = 


1 


i/3 


>V / T+F y/T+W, 

1 iy 




■ 7 2 V 1 +7 2 


(19) 


(*!*)« = [ MOMOdZ+ f MOMOdt 

Jr Jr 

and the differential operator 

Co ■ D(C 0 ) 

is self-adjoint and maximal monotoneP 9 So, it is invert¬ 
ible and we can rewrite the equation © as 

* = £o 1 [n-c I ]v. ( 20 ) 

Since D(Cq) C S and S is compactly embedded in 
Cq 1 is a compact integral operator. 


B. Exact formulae 

We assume that, for each A € R 3 , we can choose \Ea £ 
E such that the map A 1-7 \Ea is a differentiable manifold 
in H. Then, we have 
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Mutatis-mutandis , we have 

and with the same arguments, we obtain 

Therefore, 

VE m i„(A) = QllV-iAllt iII^aIH ^u^aII^ 

and 

f 1 d 

Emm (A) — E m ^ n ( 0 ) + J VE m j„(A(s)) • —A (s) ds , 

for any smooth path A(s) in K 3 joining the points (0, 0, 0) 
and (Ai, A 2 , A 3 ). In particular, for the linear path A(s) = 
sA = (sAi, SA 2 , SA 3 ), 0 < s < 1, for which we have the 
following formula 

Emin (A) = E m i ra ( 0 , 0 , 0 ) + — J ^|| V’slA II 4 A 1 

+ ||V’s 2 a||!A 2 + 2 ||' 0 s ia^ S 2 a|||A 3 ') ds . 


( 21 ) 

The chemical potentials /zi and p .2 as function of the 
parameter A can be easily calculated by multiplying the 
first equation of ([b]) by ipi\, the second by ip 2 \ and taking 
the integral over R. By this calculation, we get 
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III. VARIATIONAL APPROACH 

We consider the following trial functions: 

/o_ \ V 4 

7=7) = 7^4 i ~) exp (~T fc £ 2 ) , k = 1,2 . (24) 

By calculating the energy E with these functions, we get: 
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where, to simplify the notation, we introduced ai = 1 . 
So, by denoting /(ti,t 2 ) = E^iq,^), it is easy to see 
that /('r 1 ,r 2 ) is bounded by bellow. Indeed, if A 3 > 0, 
we have 
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and the conclusion is evident. Otherwise, notice that 
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Hence, f(T~i,T 2 ) riches its minimum at some 
Tfc(Ai, A 2 , A 3 ), (A; = 1,2) which are necessarily solu¬ 
tions of the algebraic system (i ^ j): 
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These are the equations to be solved in order to obtain 
ti(A) and 72 (A) which will be used in the formulas of 
l^app, i(A) and n a pp, 2 (A) (see below). Notice that if Ai 7 ^ 
A 2 , the respective roots are different even in the case Vq = 
0, Ni = N 2 = N and mj = m 2 - Indeed, by subtracting 
the first equation from the second one in (251, we obtain: 
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and we see that, if ti = t 2 , then Ait^ 2 — A 2 t^ ~ = 0 , 
which implies that Ai = A 2 . 


By choosing, eq = i = 1,2, the equations (25) 
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A. Approximate formulae 


Let cti(A) and 172(A) with A = (Ai, A 2 , Ai 2 , A 2 i), the 
solution of the system ([26]). 

Using Eqs. (22) and (|23[) , a direct calculation gives: 
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B. Properties of the wavefunction and the minimal 
energy 


As it was pointed out in Eq. (16), each component of 
\l/(£) in Eq. ([9]) behaves as a Gaussian as £ —► ±00, for all 
values of A, /i and Vo- In a general way, this behaviour 
justify the selection of the trial function (24). Never¬ 
theless, as it is achieved in Fig. (JT|) , the variation of the 
wavefunction of one specie with respect to the optical lat¬ 
tice intensity, Vo = Vl/Sw and the reduced wavelength, 
a -1 = d/(2ln) cannot be accounted by a Gaussian trial 


function (24). The strong variation of the optical lattice 
potential T7(£) = —Vo cos 2 (a£) with respect to a and Vo, 
keep off the contribution of the monotonic behavior of 
the harmonic potential £ 2 to order parameter. Thus, the 
variational approach presented here does not allow good 
results in the case Vo 0 is large enough. Indeed, by an 
effective numerical solution of the ID Gross-Pitaevskii 
equation we obtain the order parameter as shown 
in Figure [T] On the other hand, if we consider the equiv¬ 
alent formula of (25) for the one component BEC, we 
obtairP^ 


(7 4 + 
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For A > 0 fixed, the function er(Vo) implicitly defined by 
Eq. (29) satisfies the differential equation 
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FIG. 1: (Color online) Normalized density probability for the 
order parameter of one specie, |V , i(0 | 2 (A 3 = 0 ), for Ai = 2 , 
d/l = 0.4 and values of the laser intensity Vo = 0, 50, and 100. 
Solution of ipi (£) taken from Ref. 1241 


The choice of a test function that takes into account 
the variation shown in the figure will be treated in a 
future publication. 

Also, the presence of two-species introduces an ef¬ 
fective interaction of the unlike particle, which is de¬ 
scribed in our model by the coefficient A 3 . The effect 
of the A3 \ipj\ 2 i/>i term on the condensates is to attract 
(A 3 < 0) or to repel (A 3 > 0) the cloud probability den¬ 
sities |V ; i(£)| 2 • Thus, in the case we are dealing with a 
strong repulsive interaction, the maximum of the density 
probability lies at £ 7^ 0. Notice that the nature of our 
trial functions does not take into account the present pe¬ 
culiarity of two-species BEC. In Sec. V below we present 
a brief discussion of this effect. 


IV. PERTURBATION THEORY 


Following the result of Eq. (20), we can write the sys¬ 
tem of coupled integral equations 


/ OO 

G(£,o[M-£/]*(m', (si) 

-00 

where the kernel 

G«,f) = ( Gl( “' ) G2( “ { , ) ) , (32) 

is the solution of the differential equations £ 0 G(£, £') = 
IS(i ; — £') and I the identity matrix. In the spectral 
representation we have the Green functiorPS 
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which shows that it is increasing and blows up for a cer¬ 
tain Vq large enough. 


,* = 1,2. (33) 
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with li = A fai and <p n (z) is the harmon ic o scillator wave- 
function . 31 Thus, inserting G(£, £') in (31) we get 
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where the vectors C = (C\, C 2 ,....) and D = ( Di,D 2 ,....) 
are given by 
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and D must fulfill the non- 


To satisfy Eqs. ([34]) and (35), the vector coefficients C 

inear system of equations 
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where l r = y/ai/a 2 , ML = (n + 1/2 - /i, ; ) S nm , T and 
P(a) are matrices given elsewhere^ and S(^) is defined 
in the Appendix A. 

The above system is an infinite generalized eigenvalue 
problem for (* = 1,2) and the vector coefficients C 
and D. An efficient algorithm for solving Eqs. (36l-([37|) 
is presented in Refill. Nevertheless, it is very useful to 
carry with explicit expressions for /x* and tpi in terms of 
the leading parameters A and Vo • Assuming that the con¬ 
tribution of the non-linear terms and the optical poten¬ 
tial appearing in the system (36)-(37) are small enough 
in comparison with that of the harmonic potentials, al¬ 
lows that the vector solutions /i , C and D can be sought 
as Taylor polynomials of the parameters A and Vo- Up 
to second order terms, and solving simultaneously the 
system (36)-(371, it is possible to show that the chemical 


potentials is given by 
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Functions f(z ), g(z) and h(z) and ch(z ) are defined in 
Appendix B. 

Finally, the dimensionless order parameter, i/q, consid¬ 
ering corrections up to the first order in Ai, Ai 2 , and Vo, 
can be expressed as 

V’ P er,i M 0 + ^ j yfn2 m (rnl)2m [2-/2 
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12 
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(—l) m 2 m— 

V ( 2 m!) 


x/T+7 2 V 1 + lr 
x (a 2 ai) m exp (-a 2 ai)| ip 2m (0 ■ (39) 

, can be summed ob- 


The series, appearing in Eq. 
taining the compact solution 

Ipper, 1 Vo (0+^G(0 V2) + VoF(€,a) 

+A 12 g(e, v / i TJ 2 ) , (40) 

where T(x\ 7 ) is reported in Refill and Q(x\ z) is defined 
in the Appendix B. For the chemical potential, n per , 2 , 
and the order parameter for the second species, ip per , 2 , 
we obtain similar expressions by just changing 1 77 2 
and l r 77 1 H r in Eqs. (38) - (40). 


V. DISCUSSION OF THE RESULTS AND 
CONCLUSIONS. 

In the following we present our results and discuss the 
reliability of the two implemented methods of solution. It 
will be useful to compare the obtained analytical expres¬ 
sions with direct numerical calculations. This compari¬ 
son allows to find ranges of values of the parameters Ai, 
A 2 , A 12 A 21 and Vq where the variational approach and 
perturbation method can be implemented for the descrip¬ 
tion and predictions of the properties of the cigar-shape 
ID two-species Bose-Einstein condensates. For the nu¬ 
merical evaluation of the system ( | 20 [ ) we choose a finite 
difference method described in Ref ill 
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Ai 

FIG. 2: (Color online) Dependence of the reduced chem¬ 
ical potential g per ,i = vi/tk/j on the dimensionless self¬ 
interaction parameter Ai for the inter-particle term A 12 = 
±0.5, ±1.0, ±1.5 and A 21 = ±0.5. Values of Vo = 0, l r w 1 
and ±2 = 1 are fixed. Dashed and solid lines represent the 
analytical results from Eqs. \27\ and ( |38| >, respectively. Sym¬ 
bols correspond to the numerical solution of Eq. §. For sake 
of comparison, the limit of one component (A 3 = 0) using 
Eq. (|27[) is shown. 



Ai 

FIG. 3: (Color online) Dimensionless chemical potentials gi 
and /T 2 as a function of Ai for several species ( l r = 0.5, 1.0 
and 2.0). The same nomenclature as in Fig. © are employed. 
In the calculation we sorted Vo = 0, A 2 = 1, A 12 = 1 and 
A 21 = 0.5. 


dicated by dashed lines, while the perturbation approach, 
using Eq. ([38]), is symbolized by solid lines. Symbols rep¬ 
resent the results obtained by direct numerical evaluation 
of Eq. ©. Taking as reference the particular limit of one 
component, where A 3 = 0, as it is shown in Fig. [2] we 
observe that the influence of the inter-specie interaction 
on the chemical potential is to increase g 1 as the term 
A 12 > 0 increases, while the opposite result is achieved, 
i.e., g± decreases if Ai 2 < 0 decreases. 

The small difference seen in the figure between the per¬ 
turbation theory with respect to the variational and nu¬ 
merical solutions for Ax > 0 lies in the range of validity 
of Eq. (38). In Refl25l it is shown that the perturba¬ 
tion theory for one component reproduces quite well the 
chemical potential if |Ai| < 2. In the present case, the 
inter-species interaction plays the role as an effective non 
linear term given by Ai |-0i | 2 ± Ai 2 |^ 2 ] 2 ■ Hence, the 
range of validity of Eq. (38) as function of Ai > 0 is re¬ 
duced if Ai 2 > 0. The opposite we can argue if A i2 < 0, 
i.e., the function /iper, i(Ai) given by (38) match the vari¬ 
ational and numerical calculations in a large range of 
values of Ai > 0. Similar arguments can be performed 
for the various combination of values of the parameters 
considered in Fig. [2j 

In Fig. [3] we checked the influence of several species, 
l r = 0.5, 1, and 2, on hi and /x 2 as function of Ai with¬ 
out optical lattice, A 2 = 1 , A 2 i = 0.5, and Ai 2 = 1 . As 
might be expected, the chemical potential g 2 is almost 
constant as a function of the self-interaction term of the 
first species Ai. We note that for l r > 1 the value of the 
chemical potential hi (^ 2 ) is reduced (increased), while 
the opposite it is obtained if l r < 1. This result is ex¬ 
plained by the fact that the effective inter-species A 3 \ipi \ 2 
depends on the mass ratio l r (see Eqs. | 6 j). ([27| , (|28| and 
@). 

It can be seen that the variational approach fits very 
well the numerical calculations, but the perturbation the¬ 
ory presents some differences as Ai > 0 (Ai < 0 ) in¬ 
creases (decreases). The same argument, as it is given in 
the analysis of Fig. [2] we can argue for the dependence of 
Hper,i on Ai and l r . Nevertheless, this analysis has to be 
taken with caution. The presence of the functions f(z) 
and g(z) in Eq. (38) establishes different ranges of valid¬ 
ity for Hi(Ai) as a function of l r . Notice, that f(z) < 0 
for z > 0 } while g{z) < 0 ( g{z) > 0 ) for z > 1 (z < 1 ) 
(see Appendix B). 


A. Chemical potentials 

First, we analyze the case when the intensity of optical 
lattice is turned off, Vo = 0. Figure © shows the reduced 
chemical potential hi as a function of the dimensionless 
non-linear term Ai for the following values of the inter¬ 
species Ai 2 = ±0.5, ±1, and ±1.5. In the calculation we 
have fixed A 2 = 1, A 21 = ±0.5, and l r = 1. Variational 


B. Influence of the optical lattice 

In Fig. [4] it is shown the behavior of the chemical po¬ 
tential as function of Ai for several values of the laser 
intensity Vo, the reduced wavelength a and the A 12 pa¬ 
rameter. Solid lines represent the calculation following 
Eq. (38), dashed lines the variational approach as given 


by Eq. (27) with A i2 = 0. Symbols correspond to the 
numerical solution of Eq. (|6|) for A 12 = 0. From Fig. [4] it 


approach calculations given by Eqs. (27) and (28) are in- can be seen that Eq. (27) does not match with the per- 
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FIG. 5: (Color online) Density profile of the species BEC 
for Vo = 0. Panel (a): As a function of A 12 for Ai = 1 and 
l r = 1. Panel (b): The same as panel (a) for l r = 1.5 Panel (c): 
Varying Ai for A 12 = 1, and l r = 1.5. Panel (d): Functions 
\ipi{£)\ 2 i i=l(2) red (blue) for A 12 = 0.5 , N2/N1=.8 (solid 
line) and A 12 = 0.8, N 2 /NI -.8 (dashed lines) 3. Here Ai = 1 
and l r = 1.5 


FIG. 4: (Color online) The same as Fig. [ 2 ] for several values 
of the reduced optical lattice intensity species (Vo = 10, a = 
27r and Vo = 60,200 with a = 57r). The influence of the 
interspecies interaction is represented by solid lines. Symbols 
are the numerical solution of the Eq. (Im and dashed lines the 
variation calculation using Eq. (27). A 2 = 1 and A 21 = 0.5. 


turbation calculations neither numerical solutions. As 
Vo increases, the variational approach becomes worse, 
reflecting the choice of the trial functions (24) we have 
employed to calculate the energy. In connection with the 
perturbation theory, the agreement is satisfactory for any 
Vo less than 200 , where a small deviation from the ex¬ 
act numerical results is achieved. As it is expected, the 
influence of the unlike interspecies interaction is to in¬ 
crease the chemical potential (the opposite is obtained if 
A 12 < 0 , not shown in the figure). 


C. Miscibility of the two species 

A central issue for a description of the properties of 
multi species is the evaluation of the order parameter as a 
function of particle-particle and interspecies interaction. 
The control of the unlike particle interaction by Feshbach 
resonancJ^ allows to tune the miscibility or not of these 
structured and the challenge to create ultracold polar 
molecules. 

Figure [5] displays the spatial distribution density 


|V’per,i( 0| 2 as function of the dimensionless parameters 
A 12 (panels (a) and (b)) and Ai (panel (c)). From 
Figs.[5](a) and (b) we observe the influence of one species 
over another. The condensate is more delocalized as the 
inter-species parameter A 12 increases. Also, as the mass 
of the second species increases, the probability density 
|'i/> P ei-,i (£)| 2 spreads on the space and the maximum of 
the wavefunction decreases. The opposite is observed for 
the attractive interaction when A 12 < 0 , i.e., the density 
profile becomes more confined at £ « 0 as A i2 decreases. 
Moreover, a stronger localization occurs as the parame¬ 
ter l r increases. In other words, the system with large 
mass difference presents a more effective attraction be¬ 
tween both components, which means that it favors the 
miscibility among both species if A 12 < 0 . A compari¬ 
son between attractive and repulsive dimensionless non¬ 
linear parameter Ai is sorted in panel (c) of the figure. 
As Ai increases from 0 to 3, the density is spread is space. 
Also, for Ai large enough, the maximum of |^per,i(C)l 
is displaced by the particle-particle repulsive interaction. 
In the case of attractive interaction, Ai < 0, the maxi¬ 
mum of the order parameter i/j per ,i(0 lies at the origin. 
For sake of clarity, in panel (d) we show the influence 
of the interaction A 12 on the density profile \ip pe r,i{£)\ 
(i=1.2). Notice that the ground state is modulated by 
the repulsive interaction induced by the species 2 and 
the maximum of density probability is shifted to £ 7 ^ 0 
as A 12 increases. From the physical point of view this 
results are clear, the species 2 is expelled off the origin 
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by the first condensate. The mutual repulsion between 
the two-species affect the spatial localization of density 
profile As we stated above, this effect is driven not only 
by the values of Ai 2 , but also by the ratio of the masses 
involved in the two condensates (see Eq. (40)). 

The density distributions results of Fig. [5| indicate in 
a general way the degree of the immiscibility or phase 
separation of binary condensate due to the interspecies 
repulsion. In our case the structure is symmetric and it 
is related with the ratio of number of particles N 2 /Ni. 
These results are in complete concordance with recent 
experimental reported observations for the s 'Rb - 133 Cs 
binary condensates.^ The trial wavefunctions (24) cannot 


take into account these behaviors over the spatial distri¬ 
bution as a function of A 12 , since they are a priori located 
at the origin. 

In conclusion, we have derived simple explicit expres¬ 
sions for the chemical potentials and order parameters in 
the case of two species of non-homogeneous BEC, where 
the system is loaded in a harmonic trap potential. We 
generalize the variational method for the case of two cou¬ 
pled GP equations, showing that the obtained closed an¬ 
alytical expressions for /p (i = 1 , 2 ) represent very good 
solutions for any values of the vector A if Vq = 0. Also, 
employing the perturbation theory we are able to get 
analytical solutions for pi and the order parameter com¬ 
ponents ipi as functions of the dimensionless vector A. 
By comparison with the numerical solutions we found 
the range of validity of the Eq. (38). By the calcula¬ 
tions we show the strong dependence of p t and ^ on the 
strengths Ai, A 2 , A i2 , A 2 i and Vq. This study gives a 
very useful result establishing the universal range where 
each solution can be easily implemented. In particular, 
the dependence of the order parameter on A* and A tJ - 
(* 7 ^ j) allows to study the immiscibility of two given 
species. We should note that the variational model here 
developed can be extended to a cubic-quintic model 20 
and allows to explore the influence of quintic nonlinear 
terms on the ideal ID two coupled pure cigar-like shape 
system. 


Acknowledgements 


with H n (z) the Hermitian polynomial^-.The matrix el¬ 
ements S mn . p i(l r ) have the followings properties: 

i) l r S rn nypl (If ) — Spl;mn(^-/lr): 

ii) S 2 mO-,oo(lr), Sk 0 - 0m (lr) and Skm-oo(lr) are equal tcPl 


^2m0;00 (Zr) — 


(-1)"V(2 m)! Z 2m 
^/n(l + P r )2™m\ (1 + Z 2 ) m ’ 


(A2) 


Sk0;0 — 


(— 1 ) 2 (k + m)\ 


V^VkM.2 ^ (*±*)! (l + p r ) 


n+k+1 5 


(_l)*Sr 2^ T (Atm±l) lk+m 

TtVk'.m'. (i + I2' ) m± w ±1 

. 1 - k - m 1 + Z 2 

F{-k - to,---; , ( A 3) 


*S)cm;Oo(Zr) — 


with T (z) the Gamma functional and F(a, /?; z) the con¬ 
fluent hypergeometric function . 33 


Using the above relations it is possible to get Eqs. (38 )- 


(40). 


Appendix B: Functions 


The functions introduced in Eq. (38) are defines as: 


/(*) 

sO) 


= In 


= In 


1 1 / 2 + z 2 

2 + 2V 2(1 + z 2 ) 


1 1 / l + 2z 2 


2 V (1 +z 2 Y 


(Bl) 

(B2) 
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Appendix A: Matrix elements 

The fourth dimensional matrix S(Z r ) is defined as 

1 


mn;pl (Zr) 


tta/2 n + m + l +Pn\m\l\p\ 

/ OO 

[exp [-(1 +lr)z 2 ] H n (z ) , 

-OO 

H rn {z)Hi{l r z)H p (l r z)] dz , (Al) 


h(z) = Ei(z)—C —In 2 ; ch(z) = Chi(z)— C— In z, (B3) 

where Ei(z ) = f* c*p(x) j g ^he exponential integral, 
Chi(z) the cosine hyperbolic integral, and C the Euler’s 
constant. 


In Eq. (401 the G(£\z) is given by 


$(£;*) = 


i z t* z'-, 2 \ 

exp(-^ 2 /2) r ex P [~§2 (1 - V ) 

Z\ZttV 7T 1 / 2 J 1 — y 2 


- 1 


-dy. 

(B4) 
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